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^ ■ Abstract 

A theory is presented (and supported by numerical simulations) for phase- 
coherent reflection of light by a disordered medium which either absorbs or 
amplifies radiation. The distribution of reflection eigenvalues is shown to be 
the Laguerre ensemble of random-matrix theory. The statistical fluctuations 
of the albedo (the ratio of reflected and incident power) are computed for 
■^j- ■ arbitrary ratio of sample thickness, mean free path, and absorption or am- 

plification length. On approaching the laser threshold all moments of the 
distribution of the albedo diverge. Its modal value remains finite, however, 
' and acquires an anomalous dependence on the illuminated surface area. 
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Recent experiments on turbid laser dyes have drawn attention to the remarkable 

properties of disordered media which are optically active. The basic issue is to understand 
the interplay of phase- coherent multiple scattering and amplification (or absorption) of radi- 
ation. A quantity which measures this interplay is the albedo a, which is the power reflected 
by the medium divided by the incident power. A thick disordered slab which is optically 
passive has a — 1. Absorption leads to a < 1 and amplification to a > 1. As the amplifica- 
tion increases the laser threshold is reached, at which the average albedo becomes infinitely 
large ||. Such a generator was referred to by its inventor V. S. Letokhov as a "laser with 
incoherent feedback" |J, because the feedback of radiation is provided by random scattering 
and not by mirrors — as in a conventional laser. 

The current renewed interest in random lasers owes much to the appreciation that ran- 
domness is not the same as incoherence. Early theoretical work on this problem was based 
on the equation of radiative transfer 0, which ignores phase coherence. Zyuzin [|[] and Feng 
and Zhang |J considered interference effects on the average albedo a, averaged over different 
configurations of the scattering centra. Their prediction of a sharpening of the backscatter- 
ing peak in the angular distribution of the average reflected intensity has now been observed 
||. The other basic interference effect is the appearance of large, sample-specific fluctuations 
of the albedo around its average. These diverge faster than the average on approaching the 



laser threshold ||10|| , so that a is no longer characteristic for the albedo of a given sample. 
In the present paper we will show that, while all moments of the distribution function P(a) 
of the albedo diverge at the laser threshold, its modal value a max remains finite. The modal 
value is the value of a at which P(a) is maximal, and hence it is the most probable value 
measured in a single experiment. The diagrammatic perturbation theory of Refs. [|8|-|T0| can 
only give the first few moments of a, and hence can not determine a max . Here we develop a 
non-perturbative random-matrix theory for the entire distribution of the reflection matrix, 
from which P(a) can be computed directly. 

We contrast the two cases of absorption and amplification. In the case of absorption, 
P(a) is a Gaussian with a width 5a smaller than the average a by a factor yN, where 
iV ~ S/X 2 ^> 1 is the number of modes associated with an illuminated area S and wavelength 
A. In the case of amplification, both 5a and a increase strongly on approaching the laser 
threshold — in a manner which we compute precisely. Below threshold, the mean and modal 
value of a coincide. Above threshold, the mean is infinite while the modal value is found to 
be 



Amax = 1+0.8 7 AT. (1) 

Here 7 denotes the amplification per mean free path, assumed to be in the range iV~ 2 <C 7 <C 
1. The existence of a finite a max is due to the finiteness of the number of modes N in a surface 
area S (ignored in radiative transfer theory). Since a max scales with iV and hence with S, 
and the incident power scales with S, it follows that the reflected power scales quadratically 
rather than linearly with the illuminated area. We suggest the name "superreflection" for 
this phenomenon. To measure the albedo in the unstable regime above the laser threshold 
we propose a time-resolved experiment, consisting of illumination by a short intense pulse 
to pump the medium beyond threshold, rapidly followed by a weak pulse to measure the 
reflected intensity before spontaneous emission has caused substantial relaxation. 
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Our work on this problem was motivated by a recent paper by Pradhan and Kumar JTT 
on the case N = 1 of a single-mode waveguide. We discovered the anomalous scaling with 
area in an attempt to incorporate the effects of mode-coupling into their approach. 

We consider the reflection of a monochromatic plane wave (frequency u, wavelength A) 
by a slab (thickness L, area S) consisting of a disordered medium (mean free path I) which 
either amplifies or absorbs the radiation. We denote by a the amplification per unit length, 
a negative value of a indicating absorption. The parameter 7 = o~l is the amplification (or 
absorption) per mean free path. We treat the case of a scalar wave amplitude, and leave 
polarization effects for future study. A discrete number N of scattering channels is defined 
by imbedding the slab in an optically passive waveguide without disorder (see Fig. ^ inset). 
The number N is the number of modes which can propagate in the waveguide at frequency 
uo. The N x N reflection matrix r contains the amplitudes r mn of waves reflected into mode 
m from an incident mode n. (The basis states of r are normalized such that each carries 
unit power.) The reflection eigenvalues R n (n = 1, 2, . . . N) are the eigenvalues of the matrix 
product rr^ . The matrix r is determined by the i? n 's and by a unitary matrix U, 

Tmn X/ k Umk U n k \J Rk ■ (2) 

Note that r mn = r nm because of time-reversal symmetry. From r one can compute the 
albedo a of the slab, which is the ratio of the reflected and incident power: 

a = J2m\ r mn\ = J2kU n kU nk Rk- (3) 

For a statistical description we consider an ensemble of slabs with different configurations 
of scatterers. As in earlier work on optically passive media [I2|, we make the isotropy 



assumption that U is uniformly distributed in the unitary group. This assumption breaks 
down if the transverse dimension W of the slab is much greater than its thickness L, but is 
expected to be reasonable if W ^ L. As a consequence of isotropy, a becomes statistically 
independent of the index n of the incident mode. We further assume that the wavelength 
A is much smaller than both the mean free path I and the amplification length 1/a. The 
evolution of the reflection eigenvalues with increasing L can then be described by a Brownian 
motion process. To describe this evolution it is convenient to use the parameterization 

Rn = 1 + /i„ € (-oo,-l) U (0,oo). (4) 

The L-dependence of the distribution P(/ii, /i2, • • • //at) of the /i's is governed by the Fokker- 
Planck equation 

with initial condition lim^o-P = NYliS(fii + 1). In the single-channel case (N = 1), the 
term Y^j^i is absent and Eq. @ reduces to the differential equation studied by Pradhan and 
Kumar [fTT] , p~3|] . The multi-channel case is essentially different due to the coupling of the 
eigenvalues by the term J^j^A^j ~ /-O -1 - This term induces a repulsion of closely separated 
eigenvalues. Equation (H) with 7 = is known as the Dorokhov-Mello-Pereyra-Kumar 
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FIG. 1. Comparison between theory and simulation of the average albedo a (upper curves, 
squares) and the variance Vara (lower curves, triangles) for L/l = 1.92 (dashed curves, open 
markers) and L/l = 9.58 (solid curves, filled markers). Negative 7 corresponds to absorption, 
positive 7 to amplification. The curves are the theoretical result (0). The data points are a 
numerical simulation of a two-dimensional lattice (L = 50d and 250d, W = 51d, N = 21), averaged 
over 100 realizations of the disorder. The inset shows schematically the system under consideration. 



(DMPK) equation [pL4|,p~5]| , and has been studied extensively in the context of electronic 
conduction We have generalized it to 7 7^ 0, by adapting the approach of Ref. [|T3| to 
a non-unitary scattering matrix. 

The average a = (a) and the variance Vara = ((a — a) 2 ) of the albedo (f|) can be 
computed by first averaging U over the unitary group and then evaluating moments of the 
.Rfc's by means of Eq. 1 



]T7] . In the limit iV — > 00 we obtain the differential equations 



d _ 
l dL a = (C 

/—Vara 
dL 



-l) 2 + 2 7 a, 

4(o-l + 7) Vara + 2N~ l a(a - l) 2 . 



(6a) 
(6b) 



Corrections are smaller by a factor |7iV 2 |~ 1//2 , which we assume to be 1. Eq. ([6a] ) for the 
average albedo is an old result of radiative transfer theory [|18|]. Eq. (|6b|) for the variance 
is new. It describes the sample-specific fluctuations of the albedo due to interference of 
multiply scattered waves. Integration of Eq. (B) yields 



a = 1 - 7 + (27 -7 2 ) 1/2 tant, 
Var a = {8N cos 4 1)" 1 (47(1 - 2 7 )L/Z + 2 7 (1 + 7) - 4 7 2 cos 2t + 2 7 (1 - 7) cos At 



+ (2 - 7)" : (2 7 - 7 2 ) 1/2 [4 7 (1 - 7) sin2t 



1 -47 + 27 2 )sin4t] 



(7a) 



(7b) 



We have abbreviated t = (27 — 7 2 ) 1//2 L/l — arcsin(l — 7) 



4 



3 
"a. 0.6 

7 0.4 
0.2 



1.0 



0.8 - A 7= 



1 1 1 1 1 INI 1 


1 II 1 INI of 


a 1 1 urn 




□ 7=1.09 








- A 7=0.66 








7=0.22 








- v 7=0.072 






- 


7=0.036 
























1 1 1 



0.0 

0.01 0.1 



1 



10 



FIG. 2. Comparison between theory and simulation of the cumulative density of the variables 
fj, n (related to the reflection eigenvalues by R n = 1 + /U" 1 ). Curves are computed from the density 
( |TlD of the Laguerre ensemble; Data points are from the simulation (L = 500d = 19.2 Z, W = 151d, 
N = 63), for a single realization of the disorder. 



Plots of Eq. (0) as a function of 7 are shown in Fig. [E], for two values of L/l. (The data 
points are numerical simulations, discussed later.) In the case of absorption (7 < 0), the 
large-L limit 



1 - 7 - (7 2 - 2 7 ) 1/2 , 



Var a r 



2N 1 -7 



(8a) 
(8b) 



can be obtained directly from Eq. @ by equating the right-hand-side to zero. The limit ([|) 
is reached when L/l ^> (-y 2 — 2 r y)~ 1 / 2 . In the case of amplification (7 > 0), Eq. (|?p holds for 
L smaller than the critical length 



/(27 -7 2 )~ 1/2 arccos(7- 1) 



(9) 



at which a and Vara diverge. This is the laser threshold P,[18|. For 7 < the large-L limit 
of the probability distribution P(a) of the albedo is well described by a Gaussian, with mean 
and variance given by Eq. (|8]). (The tails are non-Gaussian, but carry negligible weight.) 
The modal value a max of the albedo equals a. For 7 > the large-L limit of P(a) can not 
be reconstructed from its moments, but needs to be determined directly. We will see that 
while a diverges, a max remains finite. 

The large-L limit -Poo^i, A*2> • • • A*iv) of the distribution of the /z's is obtained by equating 
to zero the expression between square brackets in Eq. (H). The result is 



Loo = Cnexp(-7(iV + 1)^) n \H ~ 



(10) 



with C a normalization constant. Eq. ( f[0|) holds for both positive and negative 7, but the 
support of Lqo depends on the sign of 7: All /z's have to be > for 7 > (amplification) and 
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< — 1 for 7 < (absorption). In what follows we take 7 > 0. The function (|T0|) is known 
in random-matrix theory as the distribution of the Laguerre ensemble fl9fl . The density 
p(/j) = {J2i3(fJ< ~ A*i)) °f the p's is a series of Laguerre polynomials, hence the name. For 
•yN 2 3> 1 one has asymptotically 



p(/i) = (iV/vr) (2 7 //i - 7 2 ) V2 , < /i < 2/7. (11) 



<2 



The square-root singularity at fi = is cut off in the exact density PD| , such that p = 7^ 
if /i < l/^N 2 . The cumulative density is plotted in Fig. |], together with the numerical 
simulations (discussed below). 

We seek the probability distribution of the albedo 



P(a) = (5(a - 1 - EkU^U^ 1 )) . (12) 

The average (■ • •) consists of the average of U over the unitary group followed by the av- 
erage of the /ifc's over the Laguerre ensemble. The averages can be done analytically for 
N~ 2 <C 7 <C 1 (in the continuum approximation [^TJ, i.e. by ignoring the discreteness of the 
eigenvalues), and numerically for any N, 7 (by Monte Carlo integration, i.e. by randomly 
sampling the Laguerre ensemble). The analytical result is an inverse Laplace transform, 

P(a) = — y_ ioo ^ exp[I S (a - l)/ 7 iV - 2/(*)][l + \f(s)} 2 , (13a) 

where / is an implicit function of the Laplace variable s: 



s-\f+ W 4 / + n~ 1/2 + 2(/ - v 4 / + /t 1 + 1 = o- ( l3b ) 



The continuum approximation ( fl3] ) is plotted in the inset of Fig. |^ (dashed curve). It is close 
to the exact numerical large- N result (solid curve). The modal value a max of the albedo is 
given by Eq. ([!]). The distribution P{a) drops off oc exp[— 2 / yN/(a — 1)] for smaller a and 
oc a~ 5 / 3 for larger a, so that all moments diverge. 

To test these predictions of random-matrix theory on a model system, we have carried out 
numerical simulations of the analogous electronic Anderson model with a complex scattering 
potential, using the recursive Green's function technique ||22]| . The disordered medium is 
modeled by a two-dimensional square lattice (lattice constant d, length L, width W). The 
relative dielectric constant e = e\ + 162 (relative to the value outside the disordered region) 
has a real part e\ which fluctuates from site to site between 1 ± 5s, and a non-fluctuating 
imaginary part e^- The multiple scattering of a scalar wave $ (wave number k = 2n/X) 
is described by discretizing the Helmholtz equation (V 2 + k 2 e)^ = 0. The mean free 
path I which enters in Eq. (|5p is obtained from the average albedo a — (1 + l/L)^ 1 without 
amplification (e 2 = 0). We choose k 2 = 1.5 d~ 2 , 5e = 1, leading to I = 26.1 d. The parameter 
a (and hence 7 = a I) is obtained from the analytical solution of the discretized Helmholtz 
equation in the absence of disorder (5e = 0). The complex longitudinal wavenumber k n of 
transverse mode n then satisfies the dispersion relation 

cos(M) + cos(mid/W) = 2- §(Axf) 2 (l + ie 2 ), (14) 
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FIG. 3. Comparison between theory and simulation of the cumulative probability distribution 
of the albedo (L = 500d = 19.2 1, 7 = 0.07). Solid curves are obtained by numerically averaging over 
the Laguerre ensemble; Data points are the results of the simulation, averaged over 100 realizations 
of the disorder. The three sets of data are for W = 25d, N = 10 (plusses), W = 51d, N = 21 
(triangles), and W = Wld, N = 42 (diamonds). The inset compares the continuum approximation 
( |l3| ) for P{a) (dashed) with the exact large- N limit of the Laguerre ensemble (solid). 

and leads to a = — 2iV~ 1 Im J2 n k n - The albedo @ is computed for normal incidence. 
Data points in Figs. 1-3 are the numerical results. The agreement with the curves from 
random-matrix theory is quite remarkable, given that there are no adjustable parameters. 

In summary, we have presented a random-matrix theory for the reflection matrix of a 
disordered medium with absorption or amplification. In the limit of a semi-infinite medium, 
the distribution of the reflection eigenvalues is that of the Laguerre ensemble. The corre- 
sponding distribution of the albedo is a Gaussian in the case of absorption. In the case of 
amplification, the distribution has diverging moments but a finite modal value. By ignoring 
spontaneous emission of radiation, we could examine the light reflected in response to an 
incident wave, separately from spontaneously generated light. In future work, we intend to 
include source terms in our description, to account for amplified spontaneous emission and 
the resulting relaxation of the unstable state above the laser threshold. We also plan to 
apply our approach to frequency-dependent fluctuations and to the case of diffusive, rather 
than plane-wave illumination. 

We thank A. Lagendijk, M. B. van der Mark, and D. S. Wiersma for helpful discussions. 
This work was supported by the Dutch Science Foundation NWO/FOM. 
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